First-principles theory of frozen-ion flexoelectricity 
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We demonstrate that the frozen-ion contribution to the flexoelectric coefficient is given solely in 
terms of the sum of third moments of the charge density distortions induced by atomic displacements, 
even for ferroelectric or piezoelectric materials. We introduce several practical supercell-based meth- 
ods for calculating these coefficients from first principles, and demonstrate them by computing the 
coefficients for C, Si, MgO, NaCl, SrTiOa, BaTiOa, and PbTiOa. Three important subtleties associ- 
ated with pseudopotentials, the treatment of surfaces, and the calculation of transverse components 
are also discussed. 

PACS numbers: 77.65.-j,77.90.-l-k 



Flexoelectricity (FxE) refers to the linear response of 
electric polarization to an applied strain gradient [1]. Be- 
cause a strain gradient breaks inversion symmetry, FxE 
is always symmetry-allowed, unlike piezoelectricity which 
arises only in noncentrosymmetric materials. The FxE 
effect is normally negligible on conventional length scales, 
but it may become very strong at the nanoscale, where 
huge strain gradients can significantly affect the func- 
tional properties of dielectric thin films, superlattices, 
and nanostructures. The possibility of large effects at 
the nanoscale with application to functional devices has 
caused a recent explosion of experimental interest in flex- 
oelectricity [2-8]. 

There have been remarkably few theoretical studies 
of FxE, the main difficulty being that strain gradients 
are inconsistent with translational symmetry. A classi- 
cal phenomenological theory focused on lattice-mediated 
contributions was proposed by Tagantsev [9, 10] and 
later applied to study FxE properties of dielectrics by 
Maranganti and Sharma [11]. A first attempt at a first- 
principles calculation of FxE is due to Hong et al. [12]. 
Recently, Resta [13] developed a first-principles theory 
of FxE that was, however, limited to simple elemental 
insulators such as Si, and was not implemented in prac- 
tice. Thus, unlike piezoelectricity, which is routinely cal- 
culated using modern first-principles methods in a ma- 
ture theoretical framework, the theory of FxE remains in 
a primitive state. 

In this Letter, we present a complete theory of the 
frozen- ion contributions to the FxE coefficient (EEC), 
which were not addressed in Refs. [9-11]. Working un- 
der mixed electric boundary conditions to be defined 
shortly, we demonstrate that the contribution of a given 
atom to the frozen-ion EEC is just proportional to the 
third moment of the change in charge density induced 
by its displacement. This is true for all insulating crys- 
tals, from elemental dielectrics to piezoelectrics and fer- 
roelectrics. Furthermore, we propose several practical 
supercell-based methods for extracting the EEC from 
ah initio calculations, show that these give consistent 
results, and discuss their relative advantages. We re- 



port the frozen-ion FECs for C, Si, MgO, NaCl, SrTiOg, 
BaTiOa, and PbTiOa, and discuss the trends that emerge 
from this data. Finally, we briefly discuss three impor- 
tant subtleties: (i) the issue of pseudopotential depen- 
dence; (ii) the question of "surface contributions" to the 
FxE; and (iii) the treatment of transverse components us- 
ing current-density response. The extension beyond the 
frozen-ion case, taking into account the internal lattice 
relaxations in response to strains and strain gradients, 
will be reported elsewhere. 

Theory. — Our approach here is essentially a generaliza- 
tion of the analysis introduced by Resta [13]. We consider 
an insulating crystal, fully relaxed at zero electric field 
E, and oriented such that one of its primitive reciprocal 
lattice vectors lies along x. We then identify one entire 
plane of atoms, corresponding to atom i in the home unit 
cell and its periodic images normal to x, and displace the 
entire plane rigidly by UQip in direction /3. This is done 
under electric boundary conditions in which the macro- 
scopic E continues to vanish away from the displaced 
plane. In general this induces a step in the macroscopic 
electrostatic potential, so that if done simultaneously to 
every A^'th plane of type i along x, it results in an aver- 
age E'x 7^ 0; instead what remains unchanged is the elec- 
tric displacement field D^. For this reason, we work at 
"mixed electric boundary conditions" (MEBC) in which 
we keep the macroscopic (i.e., supercell-averagcd) fields 
fixed to Ey = Ez = and Dx — 47rPs,a;, where Pg is the 
spontaneous polarization of the undeformed crystal. 

We define the planar-averaged change of charge density 
induced by this displacement to be 

dpinx + x) 



(1) 



where p{x) is the y-z planar average of p(r) and Ti is the 
location of atom i in the unit cell. We also define the 
moments of the induced charge redistribution via 



A I dx f,[i{x)x'', 



(2) 



where A is the cell area normal to x. Note that the ze- 
roth moment Q -J^'^'* vanishes due to charge conservation. 
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and that can be identified as the "Callen" or "lon- 

gitudinal" dynamical charge. 

By definition the frozen-ion FEC describes the P in- 
duced by a homogeneous strain gradient v, that is, 



(3) 



where / is a cell index and a is the lattice constant along 
X. In the spirit of Martin [15] and Resta [13], we ap- 
proach this state via the long-wave [q — 0) limit of a 
displacement wave uup = u^^e*''^'''"*''^"', where u^^ = up 
(independent of i) is small enough that a linear-response 
approach is appropriate. Then the charge density in- 
duced by the displacement of sublattice i is 



Uii3 ^ 
I 



^ fip{x -la-Tix). (4) 



This has Fourier components at Pip{q + G) at all G — 
27rm/a, but we focus on the G = component defined 
by Pipiq) = dx e~'^'^^ Pii3{x) and obtain 



Ptp[q) = 

a 



dx' I 



Uil3 I . {1.x) 



(5) 



where x' = x — la — Tix is used to obtain the first line and 
the series expansion of e~*^^ is used to obtain the second 
(terms of order q^ and higher have been dropped), and 
V = aA is the cell volume [14] . Restoring Ui/j = we get 
a total pp{q) = '^iPipiq), and using Poisson's equation 
in the form p{q) — —iqPx{q), this implies a polarization 
modulation 



PxAi) 



V 
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(6) 



where Q*'^ 



first term of Eq. (5) has dropped out due to the acoustic 



sum rule Qi^i'^'' = 0- 

Now we define the (unsymmetrized) strain tensor and 
gradient of the strain tensor to be, respectively, 



dup{r) 



drs 



(7) 



For the wave up{r) = w^e"^^ this implies ripx{q) = iqup 
and vpxxiq) = —q^up, with other elements such as rjpy 
vanishing. We also define the (unsymmetrized) frozen- 
ion piezoelectric and FxE coefficients to be 



dPo. 



which we interpret in the spirit of the long-wave method 
— liuiq^Q dPa{q)/d'qp-y{q) etc. Combining the 
above expressions with Eq. (6), it follows that [14] 

dxfix — TyTpQ B 1 (^) 



(8) 
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FIG. 1: (Color online) Supercell geometries. Large (red) and 
small (green) dots are two species of atoms; open dots are 
atoms before being displaced as shown by arrows. Rectangles 
indicate supercells. (a) Bulk supercell for Method A. (b) Slab 
supercell for method B; vertical lines indicate dipole correc- 
tion layers in vacuum, (c) Bulk supercell for Method C. 



px/SxX — QY h 



(10) 



Eq. (9) expresses the frozen-ion (or "purely elec- 
tronic") piezoelectric tensor in terms of induced 
quadrupoles quantified by the elements of Q'^^'^K This 
is basically the same as the result given in the classic 
paper of Martin [15], except that here all quantities are 
defined in the MEBC (fixed D^, Ey, and Ez). Similarly, 
Eq. (10) corresponds to the induced-octupole formula- 
tion derived in Resta's Ref. [13] and agrees with Eq. (22) 
therein (our Q^"^) is Resta's AQ'^'^^). Note, however, that 
Resta's derivation was limited to elemental (and there- 
fore non-polar and non-piezoelectric) crystals. Instead, 
the derivation here is general, showing that the frozen- 
ion FxE response has contributions only from the induced 
octupole term. 

First-principles calculations. — To compute the FECs 
from Eq. (10) using ab-initio methods, we need to set 
up a supercell calculation that allows us to calculate 
the fip{x) and, from these, the Q'fp'^\ under MEBC 
{ADx—Ey—Ez—0). We have designed three indepen- 
dent procedures to accomplish this, using three differ- 
ent supercell configurations. In Method A, shown in 
Fig. 1(a), a supercell is built from N repetitions of the 
bulk cell, and then two atomic layers are displaced in 
opposite directions under the usual boundary conditions 
in which the supercell-averaged E=0. Since the induced 
dipoles are equal and opposite, they compensate each 
other, AP = AD = 0, and the MEBC are satisfied. 
In Method B, shown in Fig. 1(b), the supercell contains 
a slab cut from the bulk material; one central layer is 
displaced, and there is an external dipole layer in the 
vacuum that is constantly readjusted so that E^ in the 
vacuum region does not change. Again, as long as there 
is no free charge on the surfaces, this enforces ADx=0. 
Finally, in Method C, illustrated in Fig. 1(c), the super- 
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TABLE I: First and third moments of displacement-induced 
charge density for MgO using three different methods. 







Q^'^ (e) 






(eBohr^) 




Method 


A 


B 


c 


A 


B 


C 


Mg 


0.63 


0.63 


0.63 


-8.91 


-8.79 


-8.35 


O 


-0.63 


-0.63 


-0.63 


-12.96 - 


-12.77 - 


-13.12 


Slim 


0.00 


0.00 


0.00 


-21.87 - 


-21. ■^)6 - 


-21.17 



cell is again bulk-like, but only one layer is displaced, now 
using a first-principles code capable of enforcing ADj.=i). 
In each case, the supercell size or slab thickness has to be 
chosen large enough that the induced charge disturbances 
fi[j{x) do not overlap or extend to the surface. 

The calculations have been performed within density- 
functional theory. We used the local-density approxi- 
mation [16] for C, Si, MgO, NaCI and SrTiOs, and the 
generalized gradient approximation [17] for BaTiOa and 
PbTiOg. We used SIESTA [18] for Methods A and B, 
ABINIT [19, 20] for Method C, and ELK [21] for the 
all-electron calculations to be discussed later. Supercells 
were built from 12 unit cells for the perovskites and 8 con- 
ventional cells for C, Si, MgO and NaCl in Method A and 
4 conventional cells for MgO in Methods B and C; slabs 
in B are separated by 20 A of vacuum. Atomic displace- 
ments of 0.04 Bohr were used in SIESTA and ABINIT, 
and 0.015 Bohr in ELK. 

Table I shows the first and third moments of MgO 
((5(2) =0 by symmetry) from Methods A-C using iden- 
tical norm-conserving pseudopotentials. Clearly the re- 
sults are in good agreement, confirming the consistent 
implementation of MEBC in all three approaches. Meth- 
ods A and B can be used to calculate FECs using stan- 
dard first-principles electronic-structure codes (although 
Method B requires a vacuum-dipole capability) , but they 
require larger supercells. Converged results can be ob- 
tained using smaller supercells with Method C, but only 
using a code that implements fixed- 2? electric boundary 
conditions [20]. 

Table II lists the moments and FECs for several mate- 
rials. For elemental and binary dielectrics, it shows that 
\lJ'xxxx\ decreases as ionicity increases. While the anion 
IQ'^^I increases from MgO to NaCl, the cation contribu- 
tion decreases, and cell volume effects also play an im- 
portant role. For all the ABO3 perovskite structures, 
the frozen-ion FECs are remarkably similar. The largest 
contribution comes from the A atoms, unlike the (Callen) 
dynamical charges Q^^\ for which Ti and Oi give domi- 
nant contributions. 

Rigid-ion model and pseudopotential dependence. — 
Note that the Q*^"^^ moments reported in Tables I and 
II, and hence the jJ-xxxx, are all negative. To see why, 
consider a model in which each cation or anion is repre- 
sented by a spherically symmetric charge pi{r) that dis- 
places rigidly as a unit. A brief calculation shows that 



TABLE IL Lattice constants (of conventional cell [14]; a and 
c for FE PbTiOa), first and third moments, and FECs as 
obtained using Method A. 





a 

(Bohr) 




(e) 


(eBohr^) 


(pC/m) 


c 


6.69 


C 


0.00 


-13.01 


-175.4 


Si 


10.22 


Si 


0.00 


-27.94 


-105.7 


MgO 


7.73 


Mg 


0.63 


-8.91 


-95.6 









-0.63 


-12.96 




NaCl 


10.66 


Na 


0.45 


-1.18 


-47.9 






CI 


-0.45 


-27.59 




SrTiOa 


7.31 


Sr 


0.39 


-54.81 


-144.7 






Ti 


1.20 


-16.48 








Oi 


-0.92 


-27.53 








O3 


-0.33 


-6.59 




BaTiOs 


7.52 


Ba 


0.40 


-65.16 


-141.9 






Ti 


1.11 


-13.80 








Oi 


-0.89 


-27.10 










— U.Ol 


— u. / 




PbTiOa 


7.43 


Pb 


0.44 


-59.03 


-156.0 






Ti 


0.83 


-25.56 








Oi 


-0.69 


-23.09 








O3 


-0.29 


-9.57 




PbTiOs 


7.35 


Pb 


0.51 


-57.40 


-148.9 


(FE) [22] 


7.88 


Ti 


0.76 


-28.41 








Oi 


-0.65 


-20.61 








O3 


-0.31 


-9.60 





Qi = J d^r x^-d^pi{r)) = Anjdrr^piir). The pos- 

itive nuclear charge at r=0 makes no contribution, so 

(3) 

within this model all Ql < 0. It is not surprising, then, 
that the real system shows a similar behavior. 

The above analysis also implies that the and 
hence Pxxxx, should depend on the treatment of the core 
density and the pseudopotential construction. (By con- 
trast, Q^^\ and hence e^xx, is unaffected.) For example, 
if the ion charge density is partitioned into core and va- 
lence contributions in the above rigid-ion model, both 
parts will contribute. We illustrate this in Table III by 
presenting results for MgO based on two approaches: an 
all-electron (AE) calculation, and a pseudopotential (PS) 
calculation in which only the change in valence electron 
density is used to define fifj{x), as for the results pre- 
sented in Tables I and II. We confirm that AE and PS 
results agree for the piezoelectric contributions, but find 
a significant difference for the FxE ones. 

This difference arises as follows. Suppose the cell- 
averaged electrostatic potentials cf)^^ and are ad- 
justed such that the valence-band maxima Evbm agree 
between the two bulk calculations. If the PS is of 
high quality, other features of the bandstructure, as 
well as forces etc., will show good agreement. However, 
_^ 0PS because —e<f){Y) is typically much deeper 
in the AE core region. Similarly, strain derivatives will 
also differ: dcj)^'^ / drj^x 7^ dcj)"^^ / drjxx- For a strain gra- 
dient at fixed Dx we have 47rAPa; = —AE^ = d(j)/dx = 
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TABLE III: Moments of MgO obtained from Method A using 
all-electron (AE) approach or pseudopotential without (PS) 
or with (PS+) rigid-core correction. 





AE 


(e) 
PS 


Q(3) 

AE 


(e Bohr^ 
PS 


) 

PS-I- 


Mg 


0.62 


0.63 


-14.57 


-8.91 


-13.76 





-0.62 


-0.63 


-12.38 - 


12.96 


-13.02 


Sum 


0.00 


0.00 


-26.95 - 


■21.87 


-26.80 



{d4>/dr]xx){drixx/dx) so that 



(d(/)/d77^^)/47r. Wc 



therefore expect n^^xx Mxxxx- Similar considerations 
apply to the theory of deformation potentials, which also 
depend on the moments Q'^^^ [23, 24]. 

The difference between Q^^.ae) g^^^ qCs.ps) jg unim- 
portant for some purposes, as for obtaining the spatial 
gradient of evBM induced by a strain gradient, where it 
cancels out of the final result. Otherwise, there is a sim- 
ple fix: for each atom type, we compute a "rigid core cor- 
rection" (RCC) gf'^CC) = 4Trfdrr'^[p^^{r) - pf^ir)] 
using the densities from free-atom AE and PS calcula- 
tions, and then add these q'^'^^*-^) corrections to the 
q(3.ps) values. We have done this for Mg and O, ob- 
taining q(3.RCC) ^ _4.85 and -0.06eBohr^ respectively. 
The corrected values, shown in the last column of Table 
III, are now in good agreement with the AE ones. 

Surface contributions. — We also considered calculating 
f^xxxx by constructing a slab supercell with two surfaces, 
as in Fig. 1(b), but applying layer displacements corre- 
sponding to the homogeneous strain gradient of Eq. (3) . 
Letting be the total slab (TS) dipole per unit area, 
we can define a FEC via n'^^x^ = Px/i^L, where v = v^xx 
and L is the slab thickness. However, wc find that pJxxx 
does not agree with the FEC computed using Methods 
A-C. On the other hand, if we compute the FEC from the 
slope of the electrostatic potential in the interior of the 
slab using window convolutions as in Ref. [13], we obtain 
IJ-xxxx = —Ex/im/ in good agreement with the results of 
Methods A-C. (In comparison with Method B, however, 
we found this method to be more diflScult to implement 
and slower to converge with slab thickness.) 

To explain why pJxxx f-i-xxxx, wc note that p^xxx 
contains contributions from the slab surfaces. To see this, 
write Anpx = <l>^ - (t>^ = Scpn — E^L — S(j)j,, where R 
and L are right and left surfaces, and for each surface 
S(j) = (j)^^"^ — (j), the difference between the vacuum level 
just outside and the macroscopic potential just inside the 
surface. Dividing by —AttvL, we find p^^xx = l^xxxx + 
((5(/)R — (5(/)L)/47ri/L. Now even if the two surfaces were 
identical initially, in the presence of the strain gradient v 
they exist at different strain states, Arj^x = >^L, and thus 
have different 64> values. In linear response we expect 
5(l)n — 6 (pi, = Ar]xx{dS4>/dr]xx), from which it follows that 
t^xxxx = t^xxxx + (d&4>/dr]xx)/'^i^- The second term is 
surface-specific [25] and reflects the dependence of the 



surface work function on local strain. 

Because we prefer that the FEC should be defined as 
a hulk property independent of surface termination, we 
adopt Pxxxx, and not f■lJxxx^ o^^' definition of the FEC. 
In a sense, Pxxxx and Pxxxx analogous respectively to 
the "proper" and "improper" contributions to piezoelec- 
tricity [26]. 

Transverse components. — The derivation ofEqs. (9-10) 
yielded e^/Sx and p-a^xx only for the case a = x. We can 
remove this restriction by replacing Eq. (1) by 



dJainx+x) 



dii, 



'Oi/3 



(11) 



where Ja{x) is the y-z planar average of the current den- 
sity in direction a induced by the adiabatic motion ■uoi/9 
of atomic plane i in direction {3, again tinder MEBC. 
Defining moments J^"'^"* = Aj dxVa.ipix) x"' , Eq. (6) 
for the polarization in direction a induced by motions in 
direction /3 is replaced by 



V 



-iqJ 



(l,x) 



q j{2,x) 



where j'^^'^^ = Y,i J^if ■ I* follows that 



M-aPxx — 2Y a/S ■ 



(12) 



(13) 



For the longitudinal case a=x, this result is equivalent to 
Eqs. (9-10), since continuity implies \7 ■Vif){r) = —fip{v), 



from which it follows that 

j(n,x) 



-(„ + l)j(».JJ. By 



contrast, the moments J„ ^'^ for a^x contain additional 

information about the transverse motions (e.g., j'^y'^ are 
transverse, or Born, charges). 

In principle, the 'Pa,ip{x) and their moments J^l^ are 
computable using the methods of density-functional per- 
turbation theory. While we have not implemented such a 
calculation here, Eq. (13) formally solves the problem of 
extending the present theory to the tensor elements Bapx 
and Pajixx ■ By carrying out similar calculations with dif- 
ferent crystal axes aligned along x, it should be possible 
to obtain the full tensors, although care must be taken 
to account for the modified interpretation of the MEBC 
after the crystal is rotated. 

Conclusions. — We have shown that the longitudinal 
frozen-ion FEC is proportional to the third moment of 
induced charge density under MEBC. An extension us- 
ing the second moment of the induced current density 
yields also the transverse FECs. This formulation is ex- 
act for all insulating crystals. Furthermore, three prac- 
tical methods for calculating FECs using ab initio meth- 
ods have been demonstrated by computing the frozen- ion 
FECs for several materials. Issues concerning pseudopo- 
tential dependence and surface effects have also been dis- 
cussed. Although it remains to include lattice contribu- 
tions associated with internal relaxations that can occur 
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in response to strains and strain gradients, our work rep- 
resents an important step in the direction of a full first- 
principles theory of FxE. 

This work was supported by ONR grant N00014-05-1- 
0054. Computations were done at the Center for Piezo- 
electrics by Design. 
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